﻿
#include "mesh.hpp"
#include <cstring>

using namespace std;


//对比source是否有needle子串，且子串从source开始
bool sub_str_cmp(const char* source, const char * needle) {
    int i = 0;
    while (needle[i] != '\0') {
        //若needle未结束，即出现不相等，则返回不相等
        if (needle[i] != source[i])
            return false;
        ++i;
    }
    //对比完成，返回相等
    return true;
}


//cphex
//private
//判断该单元是否为断层附近单元，在faults中记录本单元中的断层点；同时计算质心
//pnts为存储了该cphex中8个cpoint实际数据的指针，且索引需要匹配
//遍历该单元中8个顶点，每个点邻接另外3个点，构成3条边，可用向量a b c表示
//然后获取a b c两两夹角余弦值，若夹角小于阈值角，则判断该点为断层导致扭曲，即为断层点
void cphex::cal_fault(cpoint* pnts,int mode, int* m_faults) {

    //3条边
    cpoint a, b, c;

    //记录各点是否已经被记录为断层点
    int is_store[8];
    memset(is_store,0,sizeof(int)*8);

    //记录角落点以及相邻3点
    int four_faults_point[4];

    //分别记录两两向量之间余弦值
    double ab,ac,bc;
    //遍历8个顶点
    //不同顶点，其邻接不同
    for (int i = 0; i < 8; i++) {
        //根据不同顶点，将a、b、c更新为对应三个向量
        //同时记录相关4个点在单元中的位置
        switch (i)
        {
            case 0:
                a = pnts[mpoint[1]] - pnts[mpoint[0]];
                b = pnts[mpoint[3]] - pnts[mpoint[0]];
                c = pnts[mpoint[2]] - pnts[mpoint[0]];
                four_faults_point[0]=0;
                four_faults_point[1]=1;
                four_faults_point[2]=3;
                four_faults_point[3]=4;
                break;
            case 1:
                a = pnts[mpoint[0]] - pnts[mpoint[1]];
                b = pnts[mpoint[2]] - pnts[mpoint[1]];
                c = pnts[mpoint[3]] - pnts[mpoint[1]];
                four_faults_point[0]=1;
                four_faults_point[1]=0;
                four_faults_point[2]=2;
                four_faults_point[3]=5;
                break;
            case 2:
                a = pnts[mpoint[1]] - pnts[mpoint[2]];
                b = pnts[mpoint[3]] - pnts[mpoint[2]];
                c = pnts[mpoint[0]] - pnts[mpoint[2]];
                four_faults_point[0]=2;
                four_faults_point[1]=1;
                four_faults_point[2]=3;
                four_faults_point[3]=6;
                break;
            case 3:
                a = pnts[mpoint[0]] - pnts[mpoint[3]];
                b = pnts[mpoint[2]] - pnts[mpoint[3]];
                c = pnts[mpoint[1]] - pnts[mpoint[3]];
                four_faults_point[0]=3;
                four_faults_point[1]=0;
                four_faults_point[2]=2;
                four_faults_point[3]=7;
                break;
            case 4:
                a = pnts[mpoint[5]] - pnts[mpoint[4]];
                b = pnts[mpoint[7]] - pnts[mpoint[4]];
                c = pnts[mpoint[6]] - pnts[mpoint[4]];
                four_faults_point[0]=4;
                four_faults_point[1]=5;
                four_faults_point[2]=7;
                four_faults_point[3]=0;
                break;
            case 5:
                a = pnts[mpoint[6]] - pnts[mpoint[5]];
                b = pnts[mpoint[4]] - pnts[mpoint[5]];
                c = pnts[mpoint[7]] - pnts[mpoint[5]];
                four_faults_point[0]=5;
                four_faults_point[1]=6;
                four_faults_point[2]=4;
                four_faults_point[3]=1;
                break;
            case 6:
                a = pnts[mpoint[7]] - pnts[mpoint[6]];
                b = pnts[mpoint[5]] - pnts[mpoint[6]];
                c = pnts[mpoint[4]] - pnts[mpoint[6]];
                four_faults_point[0]=6;
                four_faults_point[1]=7;
                four_faults_point[2]=5;
                four_faults_point[3]=2;
                break;
            case 7:
                a = pnts[mpoint[6]] - pnts[mpoint[7]];
                b = pnts[mpoint[4]] - pnts[mpoint[7]];
                c = pnts[mpoint[5]] - pnts[mpoint[7]];
                four_faults_point[0]=7;
                four_faults_point[1]=6;
                four_faults_point[2]=4;
                four_faults_point[3]=3;
                break;
            default:
                break;
        }


        //获取夹角余弦值
        ab = get_cos_abs(a, b);
        ac = get_cos_abs(a, c);
        bc = get_cos_abs(c, b);

        //判断两两夹角是否大于阈值角
        if (ab < FAULT_MIN_COS || acos(ac) + acos(bc)> M_PI_4_FAULTS+0.01) {
            int related_pnts[3];
            //为faults分配空间
            if (has_faults == 0 && mode == 0) {
                faults = new int[8];
            }
            //根据不同向量夹角，将对应点存储到faults中
            related_pnts[0] = four_faults_point[0];
            related_pnts[1] = four_faults_point[1];
            related_pnts[2] = four_faults_point[2];
            if (mode == 0) {
                for (int i = 0; i < 3; i++) {
                    if (is_store[related_pnts[i]] == 0) {
                        faults[has_faults++] = mpoint[related_pnts[i]];
                        is_store[related_pnts[i]] = 1;
                    }
                }
            }
            else {
                if (acos(ac) + acos(bc) > M_PI_4_FAULTS+0.01)
                    m_faults[i] = 2;
                else
                    m_faults[i] = 1;
            }
        }
    }

    //更新点的断层属性
    for (int i = 0; mode == 0 && i < has_faults; i++) {
        pnts[faults[i]].is_fault = 1;
    }

    cal_center(pnts);
    cal_size(pnts);
}


//六面体三方向加密操作
void cphex::split(cphex *&sub,vector<cpoint> &pdata,int &psize,unordered_map<cpoint,int,cpoint_hash> &point_map){
    cpoint point_on_face[6];    //六个面上增加的点
    int id_point_on_face[6];    //六个面上增加的点的编号值
    cpoint point_on_edge[12];   //12条边上增加的点
    int id_point_on_edge[12];   //边上增加的点的编号值
    unordered_map<cpoint,int,cpoint_hash>::iterator it; //哈希迭代器
    for (int i = 0; i < 6; i++) {
        //面上的点为每个面的形心
        point_on_face[i] = (pdata[mpoint[split_face_stat[i][0]]] +
                            pdata[mpoint[split_face_stat[i][1]]] +
                            pdata[mpoint[split_face_stat[i][2]]] +
                            pdata[mpoint[split_face_stat[i][3]]]) / 4;
        //去重
        it = point_map.find(point_on_face[i]);
        if (it == point_map.end()) {
            if (pdata[mpoint[split_face_stat[i][0]]].is_fault == pdata[mpoint[split_face_stat[i][1]]].is_fault &&
                pdata[mpoint[split_face_stat[i][1]]].is_fault == pdata[mpoint[split_face_stat[i][2]]].is_fault &&
                pdata[mpoint[split_face_stat[i][2]]].is_fault == pdata[mpoint[split_face_stat[i][3]]].is_fault) {
                point_on_face[i].is_fault = pdata[mpoint[split_face_stat[i][0]]].is_fault;
            }
            id_point_on_face[i] = pdata.size();
            point_map.insert({ point_on_face[i], pdata.size() });
            pdata.push_back(point_on_face[i]);
        }
        else {
            id_point_on_face[i] = it->second;
        }
    }

    for (int i = 0; i < 12; i++) {
        //边上的点为每条边的中点
        point_on_edge[i] = (pdata[mpoint[split_edge_stat[i][0]]] + pdata[mpoint[split_edge_stat[i][1]]] )/2;
        //去重
        it = point_map.find(point_on_edge[i]);
        if (it == point_map.end()) {
            if (pdata[mpoint[split_edge_stat[i][0]]].is_fault == pdata[mpoint[split_edge_stat[i][1]]].is_fault)
                point_on_edge[i].is_fault = pdata[mpoint[split_edge_stat[i][0]]].is_fault;
            id_point_on_edge[i] = pdata.size();
            point_map.insert({ point_on_edge[i], pdata.size() });
            pdata.push_back(point_on_edge[i]);
        }
        else {
            id_point_on_edge[i] = it->second;
        }
    }
    
    int centerid=0;
    //体中心点
    cpoint t_center=(pdata[mpoint[0]] +pdata[mpoint[1]] +pdata[mpoint[2]] +pdata[mpoint[3]] +pdata[mpoint[4]] +pdata[mpoint[5]] +pdata[mpoint[6]] +pdata[mpoint[7]])/8;
    //去重
    it = point_map.find(t_center);
    if (it == point_map.end()) {
        centerid = pdata.size();
        point_map.insert({ t_center, pdata.size() });
        pdata.push_back(t_center);
    }
    else {
        centerid = it->second;
    }

    //构建8个六面体的拓扑
    for(int i=0;i<8;i++){
        sub[i].sub_init(stat,i,mpoint,id_point_on_face,id_point_on_edge,centerid);
        sub[i].center= (pdata[sub[i].mpoint[0]] + pdata[sub[i].mpoint[1]] + pdata[sub[i].mpoint[2]] + pdata[sub[i].mpoint[3]] + pdata[sub[i].mpoint[4]] + pdata[sub[i].mpoint[5]] + pdata[sub[i].mpoint[6]] + pdata[sub[i].mpoint[7]]) / 8;
    }
}

//作为加密后子六面体调用，用于构建各点拓扑
//pos为该子六面体对于原始六面体的位置
//root_mpoint为原始六面体的各顶点
//id_point_on_face为面上各点
//id_point_on_edge为边上各点
//centerid为体中心点
void cphex::sub_init(int new_stat,int pos,int *root_mpoint,int *id_point_on_face,int *id_point_on_edge,int centerid){
    stat=new_stat;
    //根据不同位置构建拓扑
    switch (pos)
    {
    case 0:
        mpoint[0]=root_mpoint[0];
        mpoint[1]=id_point_on_edge[0];
        mpoint[2]=id_point_on_face[0];
        mpoint[3]=id_point_on_edge[1];
        mpoint[4]=id_point_on_edge[4];
        mpoint[5]=id_point_on_face[1];
        mpoint[6]=centerid;
        mpoint[7]=id_point_on_face[2];
        break;
    case 1:
        mpoint[0] = id_point_on_edge[0];
        mpoint[1] = root_mpoint[1];
        mpoint[2] = id_point_on_edge[3];
        mpoint[3] = id_point_on_face[0];
        mpoint[4] = id_point_on_face[1];
        mpoint[5] = id_point_on_edge[7];
        mpoint[6] = id_point_on_face[4];
        mpoint[7] = centerid;
        break;
    case 2:
        mpoint[0] = id_point_on_edge[1];
        mpoint[1] = id_point_on_face[0];
        mpoint[2] = id_point_on_edge[2];
        mpoint[3] = root_mpoint[3];
        mpoint[4] = id_point_on_face[2];
        mpoint[5] = centerid;
        mpoint[6] = id_point_on_face[3];
        mpoint[7] = id_point_on_edge[5];
        break;
    case 3:
        mpoint[0] = id_point_on_face[0];
        mpoint[1] = id_point_on_edge[3];
        mpoint[2] = root_mpoint[2];
        mpoint[3] = id_point_on_edge[2];
        mpoint[4] = centerid;
        mpoint[5] = id_point_on_face[4];
        mpoint[6] = id_point_on_edge[6];
        mpoint[7] = id_point_on_face[3];
        break;
    case 4:
        mpoint[0]=id_point_on_edge[4];
        mpoint[1]=id_point_on_face[1];
        mpoint[2]=centerid;
        mpoint[3]=id_point_on_face[2];
        mpoint[4]=root_mpoint[4];
        mpoint[5]=id_point_on_edge[8];
        mpoint[6]=id_point_on_face[5];
        mpoint[7]=id_point_on_edge[9];
        break;
    case 5:
        mpoint[0] = id_point_on_face[1];
        mpoint[1] = id_point_on_edge[7];
        mpoint[2] = id_point_on_face[4];
        mpoint[3] = centerid;
        mpoint[4] = id_point_on_edge[8];
        mpoint[5] = root_mpoint[5];
        mpoint[6] = id_point_on_edge[11];
        mpoint[7] = id_point_on_face[5];

        break;
    case 6:
        mpoint[0] = id_point_on_face[2];
        mpoint[1] = centerid;
        mpoint[2] = id_point_on_face[3];
        mpoint[3] = id_point_on_edge[5];
        mpoint[4] = id_point_on_edge[9];
        mpoint[5] = id_point_on_face[5];
        mpoint[6] = id_point_on_edge[10];
        mpoint[7] = root_mpoint[7];

        break;
    case 7:
        mpoint[0] = centerid;
        mpoint[1] = id_point_on_face[4];
        mpoint[2] = id_point_on_edge[6];
        mpoint[3] = id_point_on_face[3];
        mpoint[4] = id_point_on_face[5];
        mpoint[5] = id_point_on_edge[11];
        mpoint[6] = root_mpoint[6];
        mpoint[7] = id_point_on_edge[10];
        break;
    default:
        break;
    }
}
//EDN cphex

//CPgrid
//public
//初始化角点单元，n_thread=n(n!=0)则多线程初始化
void CPgrid::init(const char* input_file, const char* base, int n_thread) {
    //根据base网格，标记需要的地质层界面
    vector<double> interface_z = mark_base(base, n_thread);
    //初始化网格
    init(input_file, n_thread);
    //
    int hid = 0;
    cpoint* pnt;
    cphex* hex;
    int cur = 0;

    //标记其他各层地质层界面上的点
    for (int kk = k-1; cur < interface_z.size() && kk >= 0; --kk) {
        hid = get_hex_id(0, 0, kk);
        pnt = pdata + hexdata[hid].mpoint[4];
        if (fabs(pnt->z - interface_z[cur]) < EPS) {
            for (int ii = 0; ii < i; ++ii) {
                for (int jj = 0; jj < j; ++jj) {
                    hex = get_hex(ii,jj,kk);
                    for (int p = 4; p < 8; ++p) {
                        pdata[hex->mpoint[p]].is_fault = cur + 11;
                    }
                }
            }
            ++cur;
        }
    }

    //标记最后一层地质层界面上的点
    for (int ii = 0; ii < i; ++ii) {
        for (int jj = 0; jj < j; ++jj) {
            hex = get_hex(ii, jj, 0);
            for (int p = 0; p < 4; ++p) {
                pdata[hex->mpoint[p]].is_fault = cur + 11;
            }
        }
    }
}

void CPgrid::init(const char* input_file, int n_thread)
{
    //存储input_file到filename
    int name_ptr = 0;
    for (name_ptr; name_ptr < 1024 && input_file[name_ptr] != '\0'; name_ptr++) {
        filename[name_ptr] = input_file[name_ptr];
    }
    filename[name_ptr] = '\0';
    //获取文件大小，存储至file_len
    ifstream fp(filename);
    if (!fp) {
        cout << "failed to open " << filename << endl;
        exit(1);
    }
    fp.seekg(0, ios::end);
    file_len = fp.tellg();
    fp.close();

    //存储线程参数
    if (n_thread < 0)
        n_thread = 0;
    threads_num = n_thread;
    all_threads = threads_num + 1;

    //读入网格文件，获取ijk尺寸、柱线坐标、zcorn位置等
    //cpgrid_info_mem(filename);
    cpgrid_info(filename);
    //检查给定线程参数
    thread_num_check();

    //若多线程初始化，则初始化多线程，并分发
    if (threads_num != 0)
    {
        threads.thread_init(threads_num);
        threads.run_thread(&CPgrid::deal_zcorn, this, filename);
    }
    deal_zcorn(filename);
    //等待子线程join，单线程安全
    threads.join();
    cpoint* temp = new cpoint[u_psize+1];
    memcpy(temp, pdata, sizeof(cpoint) * u_psize);
    delete[]pdata;
    pdata = temp;
    //释放中间数据空间
    delete[]ori_data;
}

//读入物理量文件
void CPgrid::read_properties(char* filename) {
    //将文件内容全部读入buff中
    char* buff;
    long long len = 0;
    ifstream fp(filename);
    if (!fp) {
        cout << "failed to open " << filename << endl;
        exit(1);
    }
    fp.seekg(0, ios::end);
    len = fp.tellg();
    buff = new char[len];

    fp.seekg(0, ios::beg);
    fp.read(buff, len);
    fp.close();

    //num标记物理量值是val的单元数目
    //参数文件中形如num*val
    //id是网格单元id值
    int num = 0, id = 0;
    double val = 0.0;
    //当前物理量数据是否结束
    int yend = 0;
    for (long long hi = 0; hi < len; hi++) {
        //youngmod

        if (sub_str_cmp(buff+hi,"YOUNGMOD\0")) {
            while (yend == 0) {
                while (!(buff[hi] >= '0' && buff[hi] <= '9')) {
                    if (buff[hi] == '/') {
                        yend = 1;
                        break;
                    }
                    ++hi;
                }
                if (yend)
                    break;
                num = get_a_num(buff, hi, num);

                val = get_a_num(buff, hi, val);
                for (int j = 0; j < num; j++) {
                    hexdata[id++].youngmode = val;
                }
            }
            id = 0;
            yend = 0;
        }
        //poissonr
        else if (sub_str_cmp(buff+hi,"POISSONR\0")) {
            while (yend == 0) {
                while (!(buff[hi] >= '0' && buff[hi] <= '9')) {
                    if (buff[hi] == '/') {
                        yend = 1;
                        break;
                    }
                    ++hi;
                }
                if (yend)
                    break;
                num = get_a_num(buff, hi, num);

                val = get_a_num(buff, hi, val);
                for (int j = 0; j < num; j++) {
                    hexdata[id++].poissonr = val;
                }
            }
            id = 0;
            yend = 0;
        }
        //ROCKDEN
        else if (sub_str_cmp(buff+hi,"ROCKDEN\0")) {
            while (yend == 0) {
                while (!(buff[hi] >= '0' && buff[hi] <= '9')) {
                    if (buff[hi] == '/') {
                        yend = 1;
                        break;
                    }
                    ++hi;
                }
                if (yend)
                    break;
                num = get_a_num(buff, hi, num);

                val = get_a_num(buff, hi, val);
                for (int j = 0; j < num; j++) {
                    hexdata[id++].rockden = val;
                }
            }
            id = 0;
            yend = 0;
        }
        //PORO
        else if (sub_str_cmp(buff+hi,"PORO\0")) {
            while (yend == 0) {
                while (!(buff[hi] >= '0' && buff[hi] <= '9')) {
                    if (buff[hi] == '/') {
                        yend = 1;
                        break;
                    }
                    ++hi;
                }
                if (yend)
                    break;
                num = get_a_num(buff, hi, num);

                val = get_a_num(buff, hi, val);
                for (int j = 0; j < num; j++) {
                    hexdata[id++].poro = val;
                }
            }
            id = 0;
            yend = 0;
        }
    }
    delete[]buff;
}

//通过单元的ijk各方向索引，获取其指针
cphex* CPgrid::get_hex(int x, int y, int z)
{
    if (x >= 0 && x < i &&
        y >= 0 && y < j &&
        z >= 0 && z < k)
    {
        return &hexdata[x + y * i + z * i * j];
    }
    else
    {
        ;// cout << "wrong index!";
        return nullptr;
    }
}

//通过给定id来获取左右前后上下邻接单元id
void CPgrid::get_neighbor(int id, int* neighbors) {
    //将单元的物理id转换为其在网格中的ijk坐标
    int x, y, z;
    x = id % i;
    y = id / i;
    z = id / (i * j);
    //获取目标单元左、右、前、后、上、下的物理id
    neighbors[0] = (x - 1 >= 0) ? id - 1 : -1;
    neighbors[1] = (x + 1 < i) ? id + 1 : -1;
    neighbors[2] = (y + 1 < j) ? id + i : -1;
    neighbors[3] = (y - 1 >= 0) ? id - i : -1;;
    neighbors[4] = (z + 1 < k) ? id + i * j : -1;
    neighbors[5] = (z - 1 >= 0) ? id - i * j : -1;
}

//根据id值，计算ijk方向上的x、y、z值
void CPgrid::get_hex_ijk(int id,int &x,int &y,int &z){
    x=id%i;
    y=id%(i*j)/i;
    z=id/(i*j);
}

//获取id值
int CPgrid::get_hex_id(int x,int y,int z) {
    if (x >= 0 && x < i &&
        y >= 0 && y < j &&
        z >= 0 && z < k)
    {
        return x + y * i + z * i * j;
    }
    else
    {
        return -1;
    }
}

//根据给定一系列单元id，来获取所构成包围盒中所有单元id
vector<int> CPgrid::get_neighbor_range(vector<int> &gid,int incr){
    vector<int> donor;
    int imin=INT_MAX,jmin=INT_MAX,kmin=INT_MAX;
    int imax=INT_MIN,jmax=INT_MIN,kmax=INT_MIN;
    int x = 0, y = 0, z = 0;
    //构建给定单元包围的包围盒
    for(int ii=0;ii<gid.size();ii++){
        get_hex_ijk(gid[ii],x,y,z);
        imin=min(imin,x);
        jmin=min(jmin,y);
        kmin=min(kmin,z);

        imax=max(imax,x);
        jmax=max(jmax,y);
        kmax=max(kmax,z);
    }
    //添加向外扩散层数
    imin = imin - incr >= 0 ? imin - incr : imin;
    jmin = jmin - incr >= 0 ? jmin - incr : jmin;
    kmin = kmin - incr >= 0 ? kmin - incr : kmin;
    imax = imax + incr <= i ? imax + incr : imax;
    jmax = jmax + incr <= j ? jmax + incr : jmax;
    kmax = kmax + incr <= k ? kmax + incr : kmax;

    //获取包围盒内所有单元的id值
    int ilen = imax - imin;
    int jlen = jmax - jmin;
    int klen = kmax - kmin;
    donor.resize((1 + ilen) * (1 + jlen) * (1 + klen), -1);
    int count = 0;
    int temp;
    for (int kk = kmin; kk <= kmax; kk++) {
        for (int jj = jmin; jj <= jmax; jj++) {
            for (int ii = imin; ii <= imax; ii++) {
                temp = get_hex_id(ii, jj, kk);
                if (temp != -1) {
                    donor[count++] = temp;
                }
            }
        }
    }
    return donor;
}

//准备adt建树所需数据接口
void CPgrid::pre_adt_info(int* ele_list, double* nodelist) {
    //传入指针需在外部做好管理
    if (ele_list == nullptr) {
        cout << "Empty adf info pointer\n";
        return;
    }

    for (int ii = 0; ii < hexsize; ii++) {
        //获取每个单元包围盒
        for (int jj = 0; jj < 8; jj++) {
            //每个单元各节点坐标值存储
            nodelist[ii * 24 + jj * 3 + 0] = pdata[hexdata[ii].mpoint[jj]].x;
            nodelist[ii * 24 + jj * 3 + 1] = pdata[hexdata[ii].mpoint[jj]].y;
            nodelist[ii * 24 + jj * 3 + 2] = pdata[hexdata[ii].mpoint[jj]].z;
        }
        //节点编号存储
        memcpy(ele_list + ii * 8, hexdata[ii].mpoint, sizeof(int) * 8);
    }
}

//准备adt建树所需数据接口
void CPgrid::pre_ele_bbox_info(double* ele_bbox) {
    if (ele_bbox == nullptr) {
        cout << "Empty adf info pointer\n";
        return;
    }
    double minx, miny, minz, maxx, maxy, maxz;
    for (int ii = 0; ii < hexsize; ii++) {
        //获取每个单元包围盒
        minx = DBL_MAX, miny = DBL_MAX, minz = DBL_MAX, maxx = DBL_MAX*-1.0, maxy = DBL_MAX*-1.0, maxz = DBL_MAX*-1.0;
        for (int jj = 0; jj < 8; jj++) {
            minx = min(pdata[hexdata[ii].mpoint[jj]].x, minx);
            miny = min(pdata[hexdata[ii].mpoint[jj]].y, miny);
            minz = min(pdata[hexdata[ii].mpoint[jj]].z, minz);

            maxx = max(pdata[hexdata[ii].mpoint[jj]].x, maxx);
            maxy = max(pdata[hexdata[ii].mpoint[jj]].y, maxy);
            maxz = max(pdata[hexdata[ii].mpoint[jj]].z, maxz);
        }
        //包围盒存储
        ele_bbox[ii * 6 + 0] = minx;
        ele_bbox[ii * 6 + 1] = miny;
        ele_bbox[ii * 6 + 2] = minz;
        ele_bbox[ii * 6 + 3] = maxx;
        ele_bbox[ii * 6 + 4] = maxy;
        ele_bbox[ii * 6 + 5] = maxz;
    }
}

//将CPgrid转换为vtk文件
void CPgrid::trans2vtk(const char* out) {
    //映射文件
    FileMapping fp(out);
    //vtk文件开头
    string buffer;
    buffer += "# vtk DataFile Version 2.0\n";
    buffer += "Unstructured Grid Example\n";
    buffer += "ASCII\n";
    buffer += "DATASET UNSTRUCTURED_GRID\n";
    buffer += "POINTS " + to_string(u_psize) + " double\n";
    fp.write(buffer);

    vector<int> fault_pnt;  //记录断层点
    string fault_scalar;    //point_data对应的标量值
    //point数据
    for (int ii = 0; ii < u_psize; ii++)
    {
        if (ii % 6 == 0) {
            buffer += '\n';
        }
        buffer += ' ' + to_string(pdata[ii].x) + ' ' + to_string(pdata[ii].y) + ' ' + to_string(pdata[ii].z);

        if(pdata[ii].is_fault!=0){
            fault_pnt.push_back(ii);
            fault_scalar += " "+to_string(pdata[ii].is_fault);
        } else {
            fault_scalar += " 0";
        }

        if(ii%10==0){
            fault_scalar+='\n';
        }
    }
    fp.write(buffer);

    //CELL数据
    int want = hexsize_valid;
    buffer += "\n\nCELLS " + to_string(hexsize_valid) + ' ' + to_string(hexsize_valid + hexsize_valid * 8) + '\n';
    int id = 0;
    for (int ii = 0; ii < hexsize_valid; ii++)
    {
        if (hexdata[ii].valid == 1)
        {
            buffer += "8 ";
            buffer += to_string(hexdata[ii].mpoint[0]) + ' '+
                      to_string(hexdata[ii].mpoint[1]) + ' '+
                      to_string(hexdata[ii].mpoint[2]) + ' '+
                      to_string(hexdata[ii].mpoint[3]) + ' '+
                      to_string(hexdata[ii].mpoint[4]) + ' '+
                      to_string(hexdata[ii].mpoint[5]) + ' '+
                      to_string(hexdata[ii].mpoint[6]) + ' '+
                      to_string(hexdata[ii].mpoint[7]) + '\n';
        }
    }

    fp.write(buffer);

    buffer += "\n\nCELL_TYPES " + to_string(want) + '\n';

    //输出CELL_TYPE 12的单元
    for (int ii = 0; ii < want; ii++) {
        if (ii % 10 == 0)
            buffer += '\n';
        buffer += "12 ";
    }

    fp.write(buffer);

    buffer = "\nCELL_DATA " + to_string(want);
    for(int i=0;i<4;i++){
        buffer+="\nSCALARS sca"+to_string(i)+" float 1\n";
        buffer+="LOOKUP_TABLE default";
        for(int hi=0;hi< want;hi++){
            if (hi % 10 == 0)
                buffer += '\n';
            if(i==0)
                buffer+=to_string(hexdata[hi].youngmode)+' ';
            else if(i==1)
                buffer+=to_string(hexdata[hi].poissonr)+' ';
            else if(i==2)
                buffer+=to_string(hexdata[hi].rockden)+' ';
            else if(i==3)
                buffer+=to_string(hexdata[hi].poro)+' ';
        }
        fp.write(buffer);
    }

    buffer.clear();
    //输出POINT_DATA以及SCALARS
    buffer+="\nPOINT_DATA "+to_string(u_psize)+'\n';
    buffer+="SCALARS isinterface int 1\n";
    buffer+="LOOKUP_TABLE default\n";
    fp.write(buffer);
    fp.write(fault_scalar);

    fp.close();
}

//输出赛题一要求属性文件
void CPgrid::out_property(char* out) {
    //文件映射
    FileMapping fp(out);
    cout<<"Writing property:"<<out<<'\n';
    string buffer;
    buffer += "YOUNGMOD";
    //杨氏模量
    fp.write(buffer);
    for (int ii = 0; ii < hexsize; ii++)
    {
        if (ii % 6 == 0) {
            buffer += '\n';
        }
        buffer += "1*" + to_string(hexdata[ii].youngmode) + ' ';
    }
    buffer += "\n/\n";
    fp.write(buffer);

    //泊松比
    buffer += "POISSONR";
    fp.write(buffer);
    for (int ii = 0; ii < hexsize; ii++)
    {
        if (ii % 6 == 0) {
            buffer += '\n';
        }
        buffer += "1*" + to_string(hexdata[ii].poissonr) + ' ';
    }
    buffer += "\n/\n";
    fp.write(buffer);

    fp.close();
}

//private
// 去重，构建的拓扑与要求的网格顶底面顺序一致
// 0为主线程，其他为子线程
// 在平行ij平面方向上，对各层并行处理
void CPgrid::deal_zcorn(char* filename, int tid )
{
    read_zcorn(tid);

    //当前线程处理的层数id范围
    int kid_start = 0;
    int kid_end = k;

    // 将数据划分到各线程
    // 余数部分交由小编号线程处理
    if (threads_num != 0)
    {
        int seg = k / all_threads;
        int nmod = k % all_threads;
        //小编号线程多处理余数部分层数
        if (tid < nmod)
        {
            kid_start = (seg + 1) * tid;
            kid_end = kid_start + (seg + 1);
        }
        else
        {
            kid_start = (nmod * (seg + 1)) + (tid - nmod) * seg;
            kid_end = kid_start + seg;
        }
    }

    // 网格退化成一个单元，去重点集大小/原点集大小=1；
    // 网格退化成两个个单元，去重点集大小/原点集大小=0.75；
    // 实际网格而言，去重点集大小/原点集大小<0.5；
    // 当前点所属单元的所在行编号
    int lineid_hex = 0;
    // ij射线编号
    int pij = 0;
    // 当前点所属单元编号
    int hexid = 0;
    // 当前点所属单元所在层编号且此处将之视为底面
    int hex_layer_id = 0;
    // 当前点在同层底面中的编号
    int hex_id_on_layer = 0;
    // 当前点所属单元在当前层的行编号
    int hex_line_id = 0;
    // 当前点所属单元在当前行中的编号
    int hex_id_on_line = 0;
    // 标记当前点是否在顶面
    int p_top = 0;
    // 标记当前点是当前面中的先给出还是后给出的连续两个
    int p_first_layer = 0;
    // 用于标记当前点在所属角点网格单元中的位置
    // bottom为单元的低z值面
    //      j   bottom:   2---3    top:   6---7
    //      ^             |   |           |   |
    //      |             |   |           |   |
    //      |             0---1           4---5
    //           ------->i
    int p_id_in_hex = 0;
    // 记录当前点与其他可能重合点对比结果
    // 若当前点与已有点重合，则check_pid为已有点编号
    int check_pid = -1;
    double x, y, z;
    bool mismatch;
    //对ori_data中各double值遍历，可根据其id，也就是其读入顺序编号，计算出该点所属网格单元id等相关信息
    for (int ic = kid_start * i * j * 8; ic < kid_end * i * j * 8; ic++)
    {
        mismatch = false;

        //线程处理到交界数据时，同步等待
        if (ic == kid_end * i * j * 8 - i * j * 4)
        {
            threads.thread_sync(tid);
        }
        hex_layer_id = ic / (i * j * 8);

        hex_id_on_layer = ic % (i * j * 4);
        hex_line_id = hex_id_on_layer / (i * 4);
        hex_id_on_line = ic % (i * 2) / 2;

        // 计算hex的编号值，hex无需去重，所以角点网格和有限元网格中编号一直
        hexid = hex_layer_id * i * j + hex_line_id * i + hex_id_on_line;

        // 计算当前点位于所属hex的位置
        p_top = (ic / (i * j * 4)) % 2;
        p_first_layer = hex_id_on_layer % (i * 4) / (i * 2);
        p_id_in_hex = 0;
        p_id_in_hex += p_top ? 4 : 0;
        p_id_in_hex += p_first_layer ? 2 : 0;
        p_id_in_hex += hex_id_on_layer % 2;

        // 每一层进位
        if (hex_id_on_layer == 0)
        {
            pij = 0;
            lineid_hex = 0;
        }
        // 每一行之间id的重复值减去
        else if (ic % (2 * i) == 0 && hex_line_id != lineid_hex)
        {
            pij -= (i + 1);
            ++lineid_hex;
        }
        z = ori_data[ic];

        // 根据Z值，以及对应Z方向控制线X Y坐标，计算出当前点X Y值
        x = getXY(ijline[pij * 2].x, ijline[pij * 2 + 1].x, ijline[pij * 2].z, ijline[pij * 2 + 1].z, z);
        y = getXY(ijline[pij * 2].y, ijline[pij * 2 + 1].y, ijline[pij * 2].z, ijline[pij * 2 + 1].z, z);

        if (ic % 2 == 0)
            ++pij;
        else if ((pij + 1) % (i + 1) == 0)
        {
            ++pij;
        }
        //跳过无效单元的计算
        if(hexdata[hexid].valid==0)
            continue;
        // 去重点，并存储，同时建立网格与点的拓扑关系
        // 因为角点网格中暗含的邻接关系，对于当前需要读入的该单元的当前点而言，根据不同位置，有不同情况
        check_pid = point_exist(hexid, p_id_in_hex, z, mismatch);
        // 在cphex中，设计的顶点邻接关系与vtk文件中CELL_TYPE 12一致，见下
        //      j   bottom:   3---2    top:   7---6
        //      ^             |   |           |   |
        //      |             |   |           |   |
        //      |             0---1           4---5
        //           ------->i
        // 故此处将2、3\6、7交换
        local_change(p_id_in_hex, 2, 3);
        local_change(p_id_in_hex, 6, 7);

        // 构建网格单元拓扑以及存储首次出现点
        if (check_pid != -1)
        {
            hexdata[hexid].mpoint[p_id_in_hex] = check_pid;
        }
        else
        {
            //线程安全
            //更新网格单元中点的索引，并保存该点
            threads.modify_var();
            hexdata[hexid].mpoint[p_id_in_hex] = u_psize;
            pdata[u_psize].z = z;
            pdata[u_psize].x = x;
            pdata[u_psize++].y = y;
            if (mismatch)
                pdata[u_psize - 1].is_fault = 1;
            if (p_id_in_hex == 6) {
                 hexdata[hexid].cal_size(pdata);
                max_edge_len = max(max_edge_len, hexdata[hexid].max_edge_len);
            }
            cpgrid_bbox[0].x = min(cpgrid_bbox[0].x, x);
            cpgrid_bbox[0].y = min(cpgrid_bbox[0].y, y);
            cpgrid_bbox[0].z = min(cpgrid_bbox[0].z, z);

            cpgrid_bbox[1].x = max(cpgrid_bbox[1].x, x);
            cpgrid_bbox[1].y = max(cpgrid_bbox[1].y, y);
            cpgrid_bbox[1].z = max(cpgrid_bbox[1].z, z);
            threads.modify_end();
        }

        //依次读取zcorn时，p_id_in_hex=6时，该点是网格单元中最后一个被读取到的点，此时判断该网格单元的断层属性
        if (p_id_in_hex == 6&&hexdata[hexid].valid) {
            hexdata[hexid].cal_fault(pdata);
        }
    }
    //线程同步
    threads.thread_sync(tid);
}

//检查给定线程数
void CPgrid::thread_num_check()
{
    //硬件支持最大线程
    int hc = thread::hardware_concurrency();
    //指定的线程数，1+子线程数
    int give = all_threads;
    //输出指定线程信息、网格k方向层数信息、硬件支持信息
    if (threads_num != 0) {
        cout << "Defined thread_num " << give << endl;;
        cout << filename<< ":CPgrid has " << k << " layers\n";
        cout << "This machine supports max " << hc << " threads\n";
    }
    //指定线程数最大不超过k且不超过硬件支持线程数的两倍
    if (all_threads > k || all_threads > 2 * hc) {
        if (all_threads > k) {
            threads_num = k - 1;
            all_threads = k;
        }
        if (all_threads > 2 * hc) {
            threads_num = 2 * hc - 1;
            all_threads = 2 * hc;
        }
    }
    //输出线程变更信息
    if (give != all_threads) {
        cout << "Change thread_num " << give << " to " << all_threads << endl;
    }
    //线程数大于硬件最大线程数，输出警告
    if (all_threads > hc) {
        cout << "Warning:Since the thread_num defined is bigger than supported, performance or result might be effeted\n";
    }
    //分配各线程ZCORN处理时指针
    if (threads_num != 0) {
        threads_data = new double* [all_threads];
        threads_data_valid = new int[all_threads];
        threads_data_size = new int[all_threads];
        int size = psize * 1.5 / all_threads;
        for (int i = 0; i < all_threads; i++) {
            threads_data[i] = new double[size];
            threads_data_size[i] = size;
            threads_data_valid[i] = 0;
        }
    }
}

//读入zcorn
void CPgrid::read_zcorn(int tid) {
    long long  zf_start = z_start_inchar;//ZCORN的Z
    long long  zf_end = z_end_inchar;//ACTNUM的A
    //划分ZCORN字段的所有字节数据
    if (threads_num != 0) {
        int period = zf_end - zf_start;
        int seg = period / all_threads;
        zf_start = (int)z_start_inchar + tid * seg;
        if (tid != threads_num)
            zf_end = zf_start + seg;
    }

    //记录当前线程读取了多少有效数据
    int dum = 0;
    int* valid = &dum;
    //存储有效数据的指针，单线程时即为ori_data
    double* mydata = ori_data;
    //记录为当前线程分配的有效数据存储的指针大小，单线程时，即为psize
    int r_size = psize;
    int* mysize = &r_size;
    //当前线程预期的有效数据数目，单线程时，即为psize
    int myseg = psize;

    if (threads_num != 0) {
        //多线程时，将各指针一一对应
        mydata = threads_data[tid];
        valid = threads_data_valid + tid;
        mysize = threads_data_size + tid;
        //在字节数据中，预期每20个字节能获取一个有效数据，得到预期有效数据数目
        myseg = (zf_end - zf_start) / 20;

    }
    //文件指针当前位置
    long long cur = zf_start;

    //判断除0号线程外，其他线程其字节数据段开头是否处于有效数据的中间
    //若处于有效数据中间，则跳过开头的不完整有效数据，交由上一线程处理，避免重复
    //注意正好处于有效数据开头的情况，额外判断一下需要
    if (tid != 0 && might_be_num(file_content[cur - 1])) {
        while (might_be_num(file_content[cur])) {
            ++cur;
        }
    }

    while (cur < zf_end) {
        //当读取的有效数据大于分配指针大小，则对其扩容
        if (*valid > (*mysize)) {
            (*mysize) *= 1.5;
            double* temp = new double[(*mysize)];
            memcpy(temp, mydata, (*valid) * sizeof(double));
            delete[]mydata;
            threads_data[tid] = temp;

            mydata = temp;
        }

        //读取一个有效数据
        //最后一个线程有可能出现读取double数据失败的情况，此时说明读到末尾，可直接返回
        mydata[*valid] = get_a_num(file_content, cur, mydata[*valid]);
        ++(*valid);
        while (!might_be_num(file_content[cur])) {
            ++cur;
        }
    }
    //线程同步
    threads.thread_sync(tid);

    //将各线程的数据汇合
    if (threads_num != 0) {
        int beg = 0;
        for (int i = 0; i < tid; i++) {
            beg += threads_data_valid[i];
        }
        memcpy(ori_data + beg, mydata, threads_data_valid[tid] * sizeof(double));
        delete[]mydata;
    }

    threads.thread_sync(tid);
    //释放多余内存
    if (threads_num != 0 && tid == 0) {
        delete[]threads_data;
        delete[]threads_data_size;
        delete[]threads_data_valid;
        if (file_content != file_mmap.read_content())
            delete[]file_content;
        else
            file_mmap.close();
    }

    return;
}

//快速读入小文件，一次性将文件读入内存
void CPgrid::cpgrid_info(char* filename)
{
    //将文件一次性读入内存
    file_mmap.init(filename, 1);
    file_content = file_mmap.read_content();
    file_len = file_mmap.len_of_reading_file();
    //标记字段状态
    int stat = 0;

    coordinate_offset[0] = 0.0;
    coordinate_offset[1] = 0.0;

    coordinate_trig_val[0] = 0.0;
    coordinate_trig_val[1] = 0.0;
    //处理file_content指针内容
    long long ptr = 0;
    while (ptr < file_len) {
        switch (file_content[ptr]) {
        case 'M': {
            if (sub_str_cmp(file_content + ptr, "MAPAXES\0")) {
                double local_coordinate[3][2];

                local_coordinate[0][0] = get_a_num(file_content, ptr, local_coordinate[0][0]);
                local_coordinate[0][1] = get_a_num(file_content, ptr, local_coordinate[0][0]);
                local_coordinate[1][0] = get_a_num(file_content, ptr, local_coordinate[0][0]);
                local_coordinate[1][1] = get_a_num(file_content, ptr, local_coordinate[0][0]);
                local_coordinate[2][0] = get_a_num(file_content, ptr, local_coordinate[0][0]);
                local_coordinate[2][1] = get_a_num(file_content, ptr, local_coordinate[0][0]);

                coordinate_offset[0] = local_coordinate[1][0];
                coordinate_offset[1] = local_coordinate[1][1];

                local_coordinate[0][0] -= local_coordinate[1][0];
                local_coordinate[0][1] -= local_coordinate[1][1];
                double all_x[2];
                all_x[0] = 1; all_x[1] = 0;
                coordinate_trig_val[0]=get_cos_2d_abs(local_coordinate[0], all_x);
                coordinate_trig_val[1] = sqrt(1.0 - coordinate_trig_val[0] * coordinate_trig_val[0]);
            }
            break;
        }
        //处理SPECGRID字段数据
        case 'S': {
            if (sub_str_cmp(file_content + ptr, "SPECGRID\0")) {
                stat += 1;
                i = get_a_num(file_content, ptr, i);
                j = get_a_num(file_content, ptr, i);
                k = get_a_num(file_content, ptr, i);

                ijsize = (i + 1) * (j + 1);
                psize = i * j * k * 8;
                hexsize = i * j * k;

                ijline = new cpoint[ijsize * 2];
                hexdata = new cphex[hexsize];
                pdata = new cpoint[psize];
                ori_data = new double[psize];
            }
            break;
        }
        //处理COORD字段数据
        case 'C': {
            if (sub_str_cmp(file_content + ptr, "COORD\0")) {
                stat += 10;
                int q = ijsize;
                int p = 0;
                while (q > 0)
                {
                    ijline[p * 2].x = get_a_num(file_content, ptr, ijline[p * 2].x);
                    ijline[p * 2].y = get_a_num(file_content, ptr, ijline[p * 2].y);
                    ijline[p * 2].z = get_a_num(file_content, ptr, ijline[p * 2].z);

                    ijline[p * 2 + 1].x = get_a_num(file_content, ptr, ijline[p * 2].z);
                    ijline[p * 2 + 1].y = get_a_num(file_content, ptr, ijline[p * 2].z);
                    ijline[p * 2 + 1].z = get_a_num(file_content, ptr, ijline[p * 2].z);
                    ++p;
                    --q;
                }
            }
            break;
        }
        //处理ZCORN字段起止信息
        case 'Z': {
            if (sub_str_cmp(file_content + ptr, "ZCORN\0")) {
                stat += 100;
                z_start_inchar = ptr;
                int look_seg_end = file_len;
                int look = file_len;
                while (z_end_inchar == 0) {
                    //每次搜索长度为200
                    look -= 200;
                    if (look < 0)
                        look = 0;

                    while (z_end_inchar == 0 && look < look_seg_end) {
                        if (file_content[look] == 'A' && sub_str_cmp(file_content + look, "ACTNUM")) {
                            z_end_inchar = look;
                            break;
                        }
                        look++;
                    }
                    //搜索区间终点是上一次起点后第20个字符
                    //以确保能搜索到ACTNUM关键字
                    look_seg_end -= 180;
                    look = look_seg_end;
                }
                ptr = z_end_inchar - 5;
            }
            break;
        }
        //处理ACTNUM字段
        case 'A': {
            if (sub_str_cmp(file_content + ptr, "ACTNUM\0")) {
                stat += 1000;
                int hex_sid = 0;
                int pre = 0;
                hexsize_valid = hexsize;
                int num = 0;
                int is_valid = 0;
                //循环读取有效单元信息
                while (ptr < file_len) {
                    num = get_a_num(file_content, ptr, num);
                    is_valid = get_a_num(file_content, ptr, num);
                    //更新对应单元有效属性
                    if (is_valid == 0) {
                        for (; hex_sid < pre + num && hex_sid < hexsize; hex_sid++) {
                            hexdata[hex_sid].valid = 0;
                        }
                        hexsize_valid -= num;
                    }
                    else {
                        hex_sid += num;
                    }
                    pre = hex_sid;
                    while (!might_be_num(file_content[ptr]) && ptr < file_len) {
                        ++ptr;
                    }
                }
            }
        }
        }
        ++ptr;
    }

    //判断角点网格文件是否符合要求
    if (stat != 111 && stat != 1111) {
        cout << "Wrong number of keyword in " << filename << endl;
        exit(1);
    }
}

int CPgrid::point_exist(int hexid, int pnt_stat, double z,bool &mismatch)
{
    int hi = hexid % (i * j) % i;
    int hj = hexid % (i * j) / i;
    int hk = hexid / (i * j);
    int ti, tj, tk;
    int target_hex_id = 0;
    int pid_in_target = 0;
    int pickit = 7 - pnt_stat;

    int temp;
    for (int q = 0; q < s_stat_len[pnt_stat]; q++)
    {
        ti = hi + s_stat[pnt_stat][q][0];
        tj = hj + s_stat[pnt_stat][q][1];
        tk = hk + s_stat[pnt_stat][q][2];
        if (ti >= 0 && ti < i &&
            tj >= 0 && tj < j &&
            tk >= 0 && tk < k)
        {
            target_hex_id = ti + tj * i + tk * i * j;
            temp = s_stat[pnt_stat][q][3];
            local_change(temp, 2, 3);
            local_change(temp, 6, 7);
            pid_in_target = hexdata[target_hex_id].mpoint[temp];
            if (pid_in_target != INT_MIN)
            {
                if (fabs(pdata[pid_in_target].z - z) < EPS)
                    return pid_in_target;
                else
                    mismatch = true;
            }
        }
    }
    return -1;
}


void CPgrid::init_neighbors(int hexid){
    int stat[6][3]={
        { 0, 0,-1},
        { 0,-1, 0},
        {-1, 0, 0},
        { 0, 1, 0},
        { 1, 0, 0},
        { 0, 0, 1}};

    int x,y,z;
    int tx,ty,tz;
    cphex *t;
    get_hex_ijk(hexid,x,y,z);
    for(int i=0;i<6;i++){
        tx=x+stat[i][0];
        ty=y+stat[i][1];
        tz=z+stat[i][2];

        t=get_hex(tx,ty,tz);

        if(tx>=0&&tx<i&&ty>=0&&ty<j&&tz>=0&&tz<k&&t!=nullptr&&t->valid){
            hexdata[hexid].neighbors[i]=t;
        }
    }
}

int fault_p = 0;
//标记地质层界面上的点
vector<double> CPgrid::mark_base(const char* basefile,int n_thread) {
    CPgrid base;
    base.init(basefile,n_thread);
    vector<double> interface_z;
    int hid = 0, pid = 0;
    cpoint* a;
    for (int k = base.k-1; k >= 0; k--) {
        hid = base.get_hex_id(0, 0, k);
        pid = base.hexdata[hid].mpoint[4];
        a = base.pdata + pid;
        interface_z.push_back(a->z);
    }
    interface_z.push_back(base.pdata[base.hexdata[hid].mpoint[0]].z);
    return interface_z;
}

//标记断层点
void CPgrid::mark_fault(const char* file){
    int id = 0;
    int neighbor[6];
    //遍历所有网格单元
    for (int ii = 0; ii < i; ii++) {
        for (int jj = 0; jj < j; jj++) {
            for (int kk = 0; kk < k; kk++) {
                id=get_hex_id(ii, jj, kk);
                neighbor[HEXUP] = get_hex_id(ii, jj, kk + 1);
                neighbor[HEXBOTTOM] = get_hex_id(ii, jj, kk - 1);
                neighbor[HEXFRONT] = get_hex_id(ii, jj - 1, kk);;
                neighbor[HEXLEFT] = get_hex_id(ii - 1, jj, kk);
                neighbor[HEXRIGHT] = get_hex_id(ii + 1, jj, kk);
                neighbor[HEXBACK] = get_hex_id(ii, jj + 1, kk);
                //对所有邻居遍历，判断邻接面上各点是否为对应同一点
                for (int hi = 0; hi < 6; hi++) {
                    if (neighbor[hi] != -1) {
                        topo_hex a(&hexdata[id]);
                        topo_hex b(&hexdata[neighbor[hi]]);
                        switch (hi) {
                        case 0:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXBOTTOM][p]) != *(b.direct[HEXUP][p])) {
                                    pdata[*(a.direct[HEXBOTTOM][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXUP][p])].is_fault = 1;
                                }
                            }
                            break;
                        case 1:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXFRONT][p]) != *(b.direct[HEXBACK][p])) {
                                    pdata[*(a.direct[HEXFRONT][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXBACK][p])].is_fault = 1;
                                }
                            }
                            break;
                        case 2:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXLEFT][p]) != *(b.direct[HEXRIGHT][p])) {
                                    pdata[*(a.direct[HEXLEFT][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXRIGHT][p])].is_fault = 1;
                                }
                            }
                            break;
                        case 3:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXBACK][p]) != *(b.direct[HEXFRONT][p])) {
                                    pdata[*(a.direct[HEXBACK][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXFRONT][p])].is_fault = 1;
                                }
                            }
                            break;
                        case 4:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXRIGHT][p]) != *(b.direct[HEXLEFT][p])) {
                                    pdata[*(a.direct[HEXRIGHT][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXLEFT][p])].is_fault = 1;
                                }
                            }
                            break;
                        case 5:
                            for (int p = 0; p < 4; p++) {
                                if (*(a.direct[HEXUP][p]) != *(b.direct[HEXBOTTOM][p])) {
                                    pdata[*(a.direct[HEXUP][p])].is_fault = 1;
                                    pdata[*(b.direct[HEXBOTTOM][p])].is_fault = 1;
                                }
                            }
                            break;
                        }
                    }
                }
            }
        }
    }

    //读取断层文件数据
    char* buff;
    long long len = 0;
    ifstream fp(file);
    if (!fp) {
        cout << "failed to open " << file << endl;
        exit(1);
    }

    fp.seekg(0, ios::end);
    len = fp.tellg();
    buff = new char[len];

    fp.seekg(0, ios::beg);
    fp.read(buff, len);
    fp.close();

    fault_bbox[0][0] = DBL_MAX;
    fault_bbox[0][1] = DBL_MAX;
    fault_bbox[0][2] = DBL_MAX;
    fault_bbox[1][0] = -1.0*DBL_MAX;
    fault_bbox[1][1] = -1.0*DBL_MAX;
    fault_bbox[1][2] = -1.0*DBL_MAX;

    double x=0.0, y=0.0, z=0.0;
    vector<cpoint> temp;
    int preid = -1;
    int curid = -1;
    vector<int> n_fault;

    for (long long ptr = 0; ptr < len;++ptr) {
        if (buff[ptr] == 'F') {
            
            //跳过第一行
            while (buff[ptr]!='\n') {
                ++ptr;
            }
            if (temp.size()) {
                fault_ij.push_back(temp);
                n_fault.push_back(fault_ij.size());
                temp.clear();
            }
            preid = -1;
            curid = -1;
        }

        if (might_be_num(buff[ptr])) {
            x = get_a_num(buff, ptr, x);
            x -= coordinate_offset[0];
            y = coordinate_offset[1] - get_a_num(buff, ptr, x);
            curid = get_a_num(buff, ptr, curid);
            z = get_a_num(buff, ptr, x);
            z *= -1.0;

            fault_bbox[0][0] = min(fault_bbox[0][0], x);
            fault_bbox[0][1] = min(fault_bbox[0][1], y);
            fault_bbox[0][2] = min(fault_bbox[0][2], z);
            fault_bbox[1][0] = max(fault_bbox[1][0], x);
            fault_bbox[1][1] = max(fault_bbox[1][1], y);
            fault_bbox[1][2] = max(fault_bbox[1][2], z);
            if (curid != preid) {
                if (preid != -1) {
                    fault_ij.push_back(temp);
                    temp.clear();
                }
                preid = curid;
            }
            temp.push_back(cpoint(x, y, z));
        }
    }

    if (temp.size()) {
        fault_ij.push_back(temp);
        n_fault.push_back(fault_ij.size());
        temp.clear();
    }
    
    
    fault_bbox[0][0] = min(fault_bbox[0][0], x);
    fault_bbox[0][1] = min(fault_bbox[0][1], y);
    fault_bbox[0][2] = min(fault_bbox[0][2], z);
    fault_bbox[1][0] = max(fault_bbox[1][0], x);
    fault_bbox[1][1] = max(fault_bbox[1][1], y);
    fault_bbox[1][2] = max(fault_bbox[1][2], z);
    delete[]buff;

    //判断各点是否位于断层文件给定数据构成的包围盒中
    for (int p = 0; p < u_psize; p++) {
        if (pdata[p].is_fault > 0 && pdata[p].is_fault <= 10) {
            if (n_fault.size() > 1) {
                cpoint f_dir;
                cpoint t_dir;
                double angle0, angle1;
                int fid = n_fault[0];
                bool checked = false;

                for (int n = 1; !checked && n < n_fault.size(); n++) {
                    for (; !checked && fid < n_fault[n] - 1; ++fid) {
                        for (int q = 0; q < fault_ij[fid].size(); q++) {
                            f_dir = fault_ij[fid + 1][q] - fault_ij[fid][q];
                            t_dir = pdata[p] - fault_ij[fid][q];
                            f_dir.z = t_dir.z;
                            angle0 = get_cos_abs(f_dir, t_dir);
                            f_dir = fault_ij[fid][q] - fault_ij[fid + 1][q];
                            t_dir = pdata[p] - fault_ij[fid + 1][q];
                            f_dir.z = t_dir.z;
                            angle1 = get_cos_abs(f_dir, t_dir);

                            if (angle0 > cos(3 / 180.0 * M_PI_4_FAULTS) || angle1 > cos(3 / 180.0 * M_PI_4_FAULTS)) {
                                pdata[p].is_fault = n + 1;
                                checked = true;
                            }
                        }
                    }
                    ++fid;
                }
            }
            if (!check_fault(pdata[p]))
                pdata[p].is_fault = 0;

        }
    }

}

bool  CPgrid::check_fault(cpoint& pnt) {
    double tol = 0;
    if (pnt.x >= fault_bbox[0][0]-tol && pnt.y >= fault_bbox[0][1]-tol && pnt.z >= fault_bbox[0][2]-tol &&
        pnt.x <= fault_bbox[1][0]+tol && pnt.y <= fault_bbox[1][1]+tol && pnt.z <= fault_bbox[1][2]+tol) {
        return true;
    }
    return false;
}

void CPgrid::add_cpoint(cpoint* t,int num) {
    if (u_psize + num > psize) {
        while(u_psize + num > psize)
            psize *= 1.5;
        cpoint* q = new cpoint[psize];
        memcpy(q, pdata, sizeof(cpoint) * u_psize);
        memcpy(q + u_psize, t, sizeof(cpoint) * num);
        delete[]pdata;
        pdata = q;
        u_psize += num;
    }
    else {
        memcpy(pdata + u_psize, t, sizeof(cpoint) * num);
        u_psize += num;
    }
}

//EDN CPgrid
void cp2msh(CPgrid& test, char* out)
{
    FileMapping fout(out);
    //ostringstream fp;
    string fp;
    cout << "writing" << endl;
    // 输出msh文件
    // gmsh的msh文件，编号从1开始
    fp = "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";
    fp+="$Nodes\n" + to_string(test.u_psize) + "\n";

    vector<string> fpnt;
    string tpnt;
    for (int i = 0; i < test.u_psize; i++)
    {
        if (test.pdata[i].is_fault!=0) {
            tpnt = " 15 1 ";
            if (test.pdata[i].is_fault == 1) {
                tpnt += " 2";
            }
            else if (test.pdata[i].is_fault == 11) {
                tpnt += " 3";
            }
            else if (test.pdata[i].is_fault == 12) {
                tpnt += " 4";
            }
            else if (test.pdata[i].is_fault == 13) {
                tpnt += " 5";
            }
            else if (test.pdata[i].is_fault == 14) {
                tpnt += " 6";
            }
            tpnt += " " + to_string(i + 1) + " ";
            fpnt.push_back(tpnt);
            //fpnt.push_back(" 15 1 "+to_string(test.pdata[i].is_fault)+" "+to_string(i));
        }
        fp+=to_string(i + 1)+" "+to_string(test.pdata[i].x)+" "+to_string(test.pdata[i].y)+" "+to_string(test.pdata[i].z)+"\n";
    }
    fout.write(fp);
    cphex *qt=test.hexdata;

//fineGRID-A
int leni=test.i;
int lenj=test.j;
int lenk=test.k;

int istart=0;
int iend=leni;
int jstart=0;
int jend=lenj;

int kstart = 0;
int kend = lenk;

int valid = 0;

    // for (int i = istart; i < iend; i++) {
    //     for(int j=jstart;j<jend;j++){
    //         for(int k=kstart;k<kend;k++){
    //             ;++valid;
    //         }
    //     }
    // }

int iistart=istart;
int iiend=iend;
int jjstart=jstart;
int jjend=jend;
int kkstart=kstart;
int kkend=kend;

    for (int i = iistart; i < iiend; i++) {
        for(int j=jjstart;j<jjend;j++){
            for(int k=kkstart;k<kkend;k++){
                ++valid;
            }
        }
    }
    fp = "$EndNodes\n$Elements\n";
    fp += to_string(valid + fpnt.size()) + '\n';

    int id = 0;
    cphex *temp;
    for (int k = kkstart; k < kkend; k++) {
        for (int i = iistart; i < iiend; i++) {
            for (int j = jjstart; j < jjend; j++) {

                temp = test.get_hex(i, j, k);
                fp += to_string(++id) + " 5 1 1 ";
                fp += to_string(temp->mpoint[0] + 1) + ' ';
                fp += to_string(temp->mpoint[1] + 1) + ' ';
                fp += to_string(temp->mpoint[2] + 1) + ' ';
                fp += to_string(temp->mpoint[3] + 1) + ' ';
                fp += to_string(temp->mpoint[4] + 1) + ' ';
                fp += to_string(temp->mpoint[5] + 1) + ' ';
                fp += to_string(temp->mpoint[6] + 1) + ' ';
                fp += to_string(temp->mpoint[7] + 1) + '\n';
            }
        }
    }

    for (string &i : fpnt) {
        fp += to_string(++id) + i + '\n';
    }

    fp +="$EndElements\n";
    // 输出msh文件END

    cout << "write end\n";
    fout.write(fp);
    fout.close();
}

void fgrid2msh(finiteGrid& fg, vector<cphex> &hexdata, vector<cpoint> &pdata, const char* out)
{
    ofstream fp;
    fp.open(out, ifstream::out);
    if (!fp) {
        cout << "Failed to open " << out << endl;
        return;
    }
    // 输出msh文件
    // gmsh的msh文件，编号从1开始
    fp << "$MeshFormat\n2.2 0 8\n$EndMeshFormat\n";
    fp << "$Nodes\n"
        << pdata.size() << "\n";

    int id = 0;
    string fpnt;
    int fc = 0;

    string str;
    for (int i = 0; i < pdata.size(); i++)
    {
        fp << i + 1 << " ";
        fp << pdata[i].x << " ";
        fp << pdata[i].y << " ";;
        fp << pdata[i].z * -1.0 << "\n";
    }

    fp<< "$EndNodes\n$Elements\n";


    int flag = 0;
    int count = 0;
    string temp;

    int is, ie, js, je, ks, ke;

    is = 0;
    ie = fg.i_withext;
    js = 0;
    je = fg.j_withext;
    ks = 0;
    ke = 1;


    int phexsize = hexdata.size();
    phexsize = 0;

    for (int i = is; i < ie; i++) {
        for (int j = js; j < je; j++) {
            for (int k = ks; k < ke; k++) {
                ++phexsize;
            }
        }
    }
    //vector<int> all=fg.children_of_ori_cphex(2, 4, 7);
    int hid = 0;
    int kkk = 1;
    for (int i = is; kkk&&i < ie; i++) {
        for (int j = js; kkk&&j < je; j++) {
            for (int k = ks; kkk&&k < ke; k++) {
                hid = fg.get_hex_id(i, j, k, fg.i_withext, fg.j_withext, fg.k_withext);
                flag = 0;
                for (int q = 0; q < 8; q++) {
                    if (hexdata[hid].mpoint[q] == INT_MIN) {
                        flag = 1;
                        ++count;
                        break;
                    }

                }

                if (flag) {
                    continue;
                }

                if (hexdata[hid].child != nullptr) {
                    for (int i = 0; i < hexdata[hid].n_child; i++) {

                        cphex* t = hexdata[hid].child + i;
                        str += to_string(++id) + " 5 1 1 ";
                        str+= to_string(t->mpoint[0] + 1) + ' ';
                        str+= to_string(t->mpoint[1] + 1) + ' ';
                        str+= to_string(t->mpoint[2] + 1) + ' ';
                        str+= to_string(t->mpoint[3] + 1) + ' ';
                        str+= to_string(t->mpoint[4] + 1) + ' ';
                        str+= to_string(t->mpoint[5] + 1) + ' ';
                        str+= to_string(t->mpoint[6] + 1) + ' ';
                        str+= to_string(t->mpoint[7] + 1) + '\n';
                    }
                }
                else {
                    str+=to_string(++id) + " 5 1 1 ";
                    str+=to_string(hexdata[hid].mpoint[0] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[1] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[2] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[3] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[4] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[5] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[6] + 1) + ' ';
                    str+=to_string(hexdata[hid].mpoint[7] + 1) + '\n';
                }
            }
        }
    }

    str = to_string(id) + "\n"+str;

    fp << str;
    fp << fpnt;
    fp << temp;

    fp << "$EndElements\n";
    // 输出msh文件END

    cout << out << " write end\n";
    fp.close();
}

void fgrid2vtk(vector<cphex> hexdata, vector<cpoint> pdata, const char* out){
    //映射文件
    FileMapping fp(out);
    //vtk文件开头
    string buffer;
    buffer += "# vtk DataFile Version 2.0\n";
    buffer += "Unstructured Grid Example\n";
    buffer += "ASCII\n";
    buffer += "DATASET UNSTRUCTURED_GRID\n";
    buffer += "POINTS " + to_string(pdata.size()) + " double\n";
    fp.write(buffer);

    vector<int> fault_pnt;  //记录断层点
    string fault_scalar = "\nPOINT_DATA " + to_string(pdata.size()) + "\nSCALARS IsInterface int 1\nLOOKUP_TABLE default\n";
    //point_data对应的标量值

    //point数据
    for (int ii = 0; ii < pdata.size(); ii++)
    {
        if (ii % 6 == 0) {
            buffer += '\n';
        }
        buffer += ' ' + to_string(pdata[ii].x) + ' ' + to_string(pdata[ii].y) + ' ' + to_string(pdata[ii].z*-1.0);

        fault_scalar += to_string(pdata[ii].is_fault)+' ';

        if (ii % 10 == 0) {
            fault_scalar += '\n';
        }
    }
    fp.write(buffer);

    //CELL数据
    int n_cell = 0;
    
    string ext_vtk;
    string ym_vtk;
    string pois_vtk;
    string rd_vtk;
    string poro_vtk;

    int id = 0;
    for (int ii = 0; ii < hexdata.size(); ii++)
    {
        if (hexdata[ii].valid == 1)
        {
            if (hexdata[ii].child != nullptr) {
                for (int i = 0; i < hexdata[ii].n_child; i++) {
                    cphex* t = hexdata[ii].child + i;
                    buffer += "8 ";
                    buffer += to_string(t->mpoint[0]) + ' ' +
                        to_string(t->mpoint[1]) + ' ' +
                        to_string(t->mpoint[2]) + ' ' +
                        to_string(t->mpoint[3]) + ' ' +
                        to_string(t->mpoint[4]) + ' ' +
                        to_string(t->mpoint[5]) + ' ' +
                        to_string(t->mpoint[6]) + ' ' +
                        to_string(t->mpoint[7]) + '\n';
                    ext_vtk += to_string(t->ext) + ' ';
                    ym_vtk += to_string(t->youngmode) + ' ';
                    pois_vtk += to_string(t->poissonr) + ' ';
                    rd_vtk += to_string(t->rockden) + ' ';;
                    poro_vtk += to_string(t->poro) + ' ';
                    ++n_cell;
                    if (n_cell % 10 == 0) {
                        ext_vtk += '\n';
                        ym_vtk += '\n';
                        pois_vtk += '\n';
                        rd_vtk += '\n';
                        poro_vtk += '\n';
                    }
                }
            }
            else {
                buffer += "8 ";
                buffer += to_string(hexdata[ii].mpoint[0]) + ' ' +
                    to_string(hexdata[ii].mpoint[1]) + ' ' +
                    to_string(hexdata[ii].mpoint[2]) + ' ' +
                    to_string(hexdata[ii].mpoint[3]) + ' ' +
                    to_string(hexdata[ii].mpoint[4]) + ' ' +
                    to_string(hexdata[ii].mpoint[5]) + ' ' +
                    to_string(hexdata[ii].mpoint[6]) + ' ' +
                    to_string(hexdata[ii].mpoint[7]) + '\n';
                ext_vtk += to_string(hexdata[ii].ext) + ' ';
                ym_vtk += to_string(hexdata[ii].youngmode) + ' ';
                pois_vtk += to_string(hexdata[ii].poissonr) + ' ';
                rd_vtk += to_string(hexdata[ii].rockden) + ' ';;
                poro_vtk += to_string(hexdata[ii].poro) + ' ';
                ++n_cell;
                if (n_cell % 10 == 0) {
                    ext_vtk += '\n';
                    ym_vtk += '\n';
                    pois_vtk += '\n';
                    rd_vtk += '\n';
                    poro_vtk += '\n';
                }
            }

        }
    }

    buffer = "\n\nCELLS " + to_string(n_cell) + ' ' + to_string(n_cell * 9) + '\n' + buffer;
    buffer += "\n\nCELL_TYPES " + to_string(n_cell) + '\n';

    //输出CELL_TYPE 12的单元
    for (int ii = 0; ii < n_cell; ii++) {
        if (ii % 10 == 0)
            buffer += '\n';
        buffer += "12 ";
    }
    fp.write(buffer);

    ext_vtk = "\nCELL_DATA " + to_string(n_cell) + "\nSCALARS CellZone int 1\nLOOKUP_TABLE default\n" + ext_vtk;
    ym_vtk = "\nSCALARS YoungMod float 1\nLOOKUP_TABLE default\n" + ym_vtk;
    pois_vtk = "\nSCALARS Poissonr float 1\nLOOKUP_TABLE default\n" + pois_vtk;
    rd_vtk = "\nSCALARS Rockden float 1\nLOOKUP_TABLE default\n" + rd_vtk;
    poro_vtk = "\nSCALARS Poro float 1\nLOOKUP_TABLE default\n" + poro_vtk;

    fp.write(ext_vtk);
    fp.write(ym_vtk);
    fp.write(pois_vtk);
    fp.write(rd_vtk);
    fp.write(poro_vtk);
    fp.write(fault_scalar);

    //输出SCALARS
    cout << out << " write end\n";
    fp.close();
}